#!/usr/bin/perl
use Getopt::Std;
#getopt("a:c:d:n:m", \%args);
getopt("a:d:n:m", \%args);
#$ambi = $args{a};
#$ambi ="chr8q24_GBM.haps.gz";
$ambi ="chr8q24_GBM.haps.gz";
chomp $ambi;
$dosage = "impute_8q24_out_ambi";
chomp $dosage;
if(!(-e $ambi))
{
	die "$ambi file not exists\n";
}
open(SAM,"gunzip -c $ambi |") or die "no file found $ambi\n";
open(DOS,"$dosage") or die "no file found $dosage\n";
#reading all the ambiguous snp positions in to hash
while($line = <SAM>)
{
	chomp $line;
	$line =~ s/\n//g;
	$line =~ s/\r//g;
	#@array_sam = split(/\s/,$line);
	@array_haps =split(/\s/,$line);
     if(($array_haps[3] eq "A" && $array_haps[4] eq "T" ) || ($array_haps[3] eq "T" && $array_haps[4] eq "A")||($array_haps[3] eq "C" && $array_haps[4] eq "G")||($array_haps[3] eq "G" && $array_haps[4] eq "C"))
     {	
	undef(@array_sam);
	push(@array_sam,"8");
	push(@array_sam,$array_haps[1]);
	push(@array_sam,"0");
	push(@array_sam,$array_haps[2]);
	$all_1=$array_haps[3];
	$all_2=$array_haps[4];
	for($z=5;$z<@array_haps;$z++)
	{
		if($array_haps[$z] eq "0")
		{
			push(@array_sam,$all_1);
		}
		else
		{
			push(@array_sam,$all_2);
		}
	}
	#print $array_sam[3]."\n";
	$line=join(" ",@array_sam);
	#die "$line\n";
	$hash_sam{$array_sam[3]} =$line;
     }	
}

undef(@array_sam);
$num_samp = 0;
$S1=0;
$S2=0;
$S3=0;
#reading the dosage file
while($line = <DOS>)
{
	$S1=0;
	$S2=0;
	$S3=0;	
	$S4=0;
	$S5=0;	
	chomp($line);
	@array = split(" ",$line);
	if(length($array[3]) == 1 && length($array[4]) == 1 && $array[4] =~ m/A|T|G|C/ && $array[3] =~ m/A|T|G|C/)
	{
		if($num_samp ==0)
		{
			#START: create header to simulate beagle dosage file
			$num_samp = @array;
			$header = "position marker alleleA alleleB";
			#print "$num_samp\n";
			$num_samp =($num_samp-5)/3;
			for($i=1;$i<=$num_samp;$i++)
			{
				$header=$header." $i $i $i";
			}
			print WRBUFF1 $header."\n";
			#STOP: create header to simulate beagle dosage file
		}
		if(exists($hash_sam{$array[2]}))
		{
			#print "$line\n";
			$line1 = $hash_sam{$array[2]};
			@array_sam = split(/\s/,$line1);
			$a1 = $array[3];
			$a2 = $array[4];
			$count_fwd = 0;
			$count_rev = 0;
			shift(@array_sam);
			shift(@array_sam);
			shift(@array_sam);
			shift(@array_sam);
			$line1 = join(" ",@array_sam);
			$j = 0;
			#print "$line\n$line1\n";
			#count the maf and checking should be greater than 30%
			$maf_check = $line1;
			
			if($line1 =~ m/A/)
			{
				$allele1_count_1 =$maf_check=~ tr/A//;
				$allele1_count_2 =$maf_check=~ tr/T//;
			#	$print = "A T";
			}
			else{
				$allele1_count_1 =$maf_check=~ tr/G//;
				$allele1_count_2 =$maf_check=~ tr/C//;
			#	$print ="G C";
			}
			#print "$print $line1 $allele1_count_1 $allele1_count_2\n";
			if($allele1_count_1 < $allele1_count_2)
			{
				$cutoff = $allele1_count_1/($allele1_count_2+$allele1_count_1);
				#print "$print $allele1_count_1 $allele1_count_2 $cutoff\n";
			}
			else
			{
				$cutoff = $allele1_count_2/($allele1_count_2+$allele1_count_1);
				#print "$print $allele1_count_2 $allele1_count_1 $cutoff\n";
			}
			#if($array[2] eq "126325828")
			#{
			#	die "test $cutoff\n";
			#}
			#if($cutoff > 0.3)
			{
				#print "counting the matches of homozygotes in both strands\n";
				for($i=5;$i<@array;$i++)
				{
					$d1 = $array[$i++];
					$d2 = $array[$i++];
					$d3 = $array[$i];
					$gnt1 = $array_sam[$j++]." ".$array_sam[$j++];
					$gnt1_flip = $gnt1;
					$gnt1_flip =~ tr/ATGC/TACG/;
					if($d1 >0.5)
					{
						$gnt2 = "$a1 $a1";
					}
					elsif($d2 > 0.5)
					{
						$gnt2 ="$a1 $a2";
					}
					elsif($d3 > 0.5)
					{
						$gnt2 = "$a2 $a2";
					}
					else
					{
						$gnt2 = "0 0";
					}
					@g1 = split(" ",$gnt1);
					@g1_f = split(" ",$gnt1_flip);
					@g2 = split(" ",$gnt2);
					#print "gnt1:$gnt1\t@g1\t\t\t\tgnt-flip:$gnt1_flip\t@g1_f\t\t\t\t\tgnt2:$gnt2\t@g2\n";
					#print "$gnt1\n";
					if($g1[0] eq $a1 && $g1[1] eq $a1)
					{
						$S1++;
					}
					elsif(($g1[0] eq $a1 && $g1[1] eq $a2 )||($g1[1] eq $a1 && $g1[0] eq $a2 ))
                                        {
                                                $S2++;
                                        }
					elsif($g1[0] eq $a2 && $g1[1] eq $a2)
                                        {
                                                $S3++;
                                        }
					else
					{
						$S4++;
					}

					if($g1[0] eq $g1[1] && $gnt1 ne "0 0" && $g2[0] eq $g2[1] && $gnt2 ne "0 0" && $gnt1 eq $gnt2)
					{
						$count_fwd++;
					}
					if($g1_f[0] eq $g1_f[1] && $gnt1_flip ne "0 0" && $g2[0] eq $g2[1] && $gnt2 ne "0 0" && $gnt1_flip eq $gnt2)
					{
						$count_rev++;
					}
				}
				undef(@array_sam);
				#fliping the input genotype according to counts
				if($count_fwd > $count_rev)
				{
					$finalgnt = $line1;
					#die  "$array[0] $array[1] $array[2] $array[3] $array[4]\n";
				}
				else
				{
					$line1 =~ tr/ATGC/TACG/;
					$finalgnt = $line1;
					#die  "$array[0] $array[1] $array[2] $array[3] $array[4]\n";
				}	
				$S5=$S1+$S2+$S3+$S4;
				print "$array[1] $array[2] $array[3] $array[4] $count_fwd $count_rev $S1 $S2 $S3 $S4 $S5\n";
				#replacing the doage with selected forward or reverse original genotypes
			}	
		}
	}
}
